Analyse spatiale et territoriale du logement social

Formation Carthageo-Geoprisme 2022

Introduction

Données RP 2018

Définir le sujet

Soit le sujet : Logements sociaux et qualification des chefs de ménages

Définir les “logements sociaux”

Logements HLM ? Logements SRU ?

Définir la notion de “qualification” ?

Le diplôme le plus élevé ? le nombre d’années d’étude ?

Définir la date

Année 2018 uniquement ? Résultats du RP 2018 (2016-2020) ?

Formuler des questions ou des hypothèses

Qu’elles soient justes ou fausses, les hypothèses permettent de cadrer l’analyse.

Diplôme et logement social

Les logements sociaux sont réservés aux ménages les moins diplômés

Âge et logement social

Les logements sociaux sont réservés aux jeunes ménages

Logement social et territoire

Les logements sociaux sont concentrés dans certains quartiers

Logement social, âge et diplômes

Les personnes diplômés quittent les logements sociaux dès que leurs revenus progressent

Organiser le travail

Sutout dans le cadre d’un groupe !

Ne collecter que les données utiles pour répondre aux questions posées

Afin de ne pas être tenté de partir dans toutes les directions

Archiver soigneusement les programmes et les résultats

Afin de pouvoir reproduire ultérieurement les analyses sur une autre période, un autre territoire

Ne pas attendre d’avoir accumulé tous les résultats pour les commenter

Car l’analyse peut suggérer des erreurs ou ouvrir de nouvelles pistes.

Partir des questions et non pas des outils

Faute de quoi on va trouver des réponses (42 …) sans savoir quelle est la question.

Charger les données statistiques

programme

tab_ind<-readRDS("data/menag2018.RDS")

résultat

FALSE   COMMUNE   ARM      IRIS ACHL AEMM
FALSE 1   75056 75101 751010101 C114 2014
FALSE 2   75056 75101 751010101  A11 2008

Préparation de l’analyse

  • Soit la relation entre logement en HLM (Y) et Diplôme le plus élevé du chef de ménage (X). Il s’agit de deux variables catégorielles (= qualitatives) que l’on va typiquement mettre en relation à l’aide d’un tableau de contingence et d’un test du chi-2. L’analyse statistique est simple sous R mais il faut tenir compte de trois difficultés

  • Le choix de la population de référence est important. Ici on va sélectionner les ménages dont la personne de référence est âgée de 25-39 ans

  • Le choix de l’espace de référence est par exemple celui du 13e arrondissement et des communes voisines d’Ivry, Gentilly et Kremlin-Bicêtre

  • la sélection ou le regroupement des diplômes est également important car cela va influer sur les résultats du test.

  • la pondération des individus doit également être prise en compte puisque le recensement est basé sur un sondage

Sélection des individus et des variables

programme

#table(tab_ind$AGEMEN8)
tab_sel<- tab_ind %>% 
  filter(AGEMEN8 == "25") %>%
  filter(COMMUNE %in% c(75113, 94037, 94041, 94043))%>%
  select(COMMUNE, DIPLM,HLML, IPONDL) 

résultats

COMMUNE DIPLM HLML IPONDL
94037 13 1 4.681816
94037 02 1 4.681816
94037 15 1 4.681816
94037 18 1 1.501089

Recodage des modalités

On cherche le code des modalités CS1 ezt HLML dans le fichier des métadonnées

meta<-readRDS("data/menag2018_meta.RDS")
metasel <- meta %>% filter(COD_VAR %in% c("DIPLM", "HLML"))
kable(metasel[,c(1,3,4)])
COD_VAR COD_MOD LIB_MOD
DIPLM 01 Pas de scolarité ou arrêt avant la fin du primaire
DIPLM 02 Aucun diplôme et scolarité interrompue à la fin du primaire ou avant la fin du collège
DIPLM 03 Aucun diplôme et scolarité jusqu’à la fin du collège ou au-delà
DIPLM 11 CEP (certificat d’études primaires)
DIPLM 12 BEPC, brevet élémentaire, brevet des collèges, DNB
DIPLM 13 CAP, BEP ou diplôme de niveau équivalent
DIPLM 14 Baccalauréat général ou technologique, brevet supérieur, capacité en droit, DAEU, ESEU
DIPLM 15 Baccalauréat professionnel, brevet professionnel, de technicien ou d’enseignement, diplôme équivalent
DIPLM 16 BTS, DUT, Deug, Deust, diplôme de la santé ou du social de niveau bac+2, diplôme équivalent
DIPLM 17 Licence, licence pro, maîtrise, diplôme équivalent de niveau bac+3 ou bac+4
DIPLM 18 Master, DEA, DESS, diplôme grande école niveau bac+5, doctorat de santé
DIPLM 19 Doctorat de recherche (hors santé)
DIPLM YY Hors résidence principale
HLML 1 Logement appartenant à un organisme HLM
HLML 2 Logement n’appartenant pas à un organisme HLM
HLML Y Hors résidence principale

s

On recode les modalités des deux variables en regroupant certaines CSP

programme

tab_sel$HLML<-as.factor(tab_sel$HLML)
levels(tab_sel$HLML)<-c("HLM-O","HLM-N",NA)
tab_sel$DIPLM<-as.factor(tab_sel$DIPLM)
levels(tab_sel$DIPLM) <- c("< BAC","< BAC","< BAC","< BAC","< BAC","< BAC",
                         "BAC","BAC",
                         "BAC+123","BAC+123","> BAC+3","> BAC+3",NA)
table(tab_sel$DIPLM)
FALSE 
FALSE   < BAC     BAC BAC+123 > BAC+3 
FALSE    1352     898    1572    2277

résultats

COMMUNE DIPLM HLML IPONDL
94037 < BAC HLM-O 4.681816
94037 < BAC HLM-O 4.681816
94037 BAC HLM-O 4.681816

Création du tableau de contingence non pondéré (FAUX)

La solution la plus simple semble être l’instruction table()

programme

tab_cont<-table(tab_sel$HLML,tab_sel$DIPLM)

résultats

< BAC BAC BAC+123 > BAC+3 Sum
HLM-O 662 410 476 307 1855
HLM-N 690 488 1096 1970 4244
Sum 1352 898 1572 2277 6099

Création du tableau de contingence pondéré (JUSTE)

On pondère avec wtd.table() du package questionr.

programme

library(questionr)
tab_cont_wtd<-wtd.table(tab_sel$HLML,tab_sel$DIPLM,
                        weights = tab_sel$IPONDL)

résultats

< BAC BAC BAC+123 > BAC+3 Sum
HLM-O 1500 928 1085 680 4192
HLM-N 1754 1220 2661 4773 10408
Sum 3254 2148 3745 5453 14600

Comparaison des niveaux de dépendance automobile

  • Tableau non pondéré … légèrement faux !
< BAC BAC BAC+123 > BAC+3 Ensemble
HLM-O 49 45.7 30.3 13.5 30.4
HLM-N 51 54.3 69.7 86.5 69.6
Total 100 100.0 100.0 100.0 100.0
  • Tableau pondéré … juste !
< BAC BAC BAC+123 > BAC+3 Ensemble
HLM-O 46.1 43.2 29 12.5 28.7
HLM-N 53.9 56.8 71 87.5 71.3
Total 100.0 100.0 100 100.0 100.0

Visualisation du tableau de contingence

On choisit l’orientation du tableau et on l’affiche avec plot()

mytable<-wtd.table(tab_sel$DIPLM,tab_sel$HLML,weights = tab_sel$IPONDL)
plot(mytable)

Visualisation améliorée du tableau de contingence

Tant qu’à faire, on améliore la figure avec des paramètres supplémentaires :

Test du Chi-deux

Ce test se réalise facilement sur le tableau de contingence avec l’instruction chisq.test() :

mytest<-chisq.test(mytable)
mytest

    Pearson's Chi-squared test

data:  mytable
X-squared = 1403.9, df = 3, p-value < 2.2e-16

Visualisation des résidus

Lorsque la relation est significative, on visualise les cases les plus exceptionnelles avec mosaicplot( …, shade = T)

Localisation aréale (sf)

Le format sf (spatial features)

La cartographie et plus généralement les opérations géométriques sur des données spatiales dans R peuvent facilement être effectuées avec le package sf (spatial features) qui crée des objets uniques rassemblant à la fois

  • un tableau de données (l’équivalent du fichier .dbf)
  • une géométrie (l’équivalent du fichier .shp)
  • une projection (l’équivalent du fichier .prj)

Lorsqu’on récupère des fonds de carte au format shapefile (.shp) ou dans d’autres formats standards comme GeoJson, la première tâche consiste donc à les convertir au formt sf afin de pouvoir les utiliser facilement dans R. L’importation se fait à l’aide de l’instruction st_read en indiquant juste le nom du fichier .shp à charger. Les autres fichiers (.dbf ou .proj) seront lus également et intégrés dans l’objet qui hérite de la double classe data.frame et sf.

Etapes de préparation des données

Dans notre exemple, nous allons suivre les étapes suivantes :

  1. Préparer les données statistiques par IRIS dans un data.frame
  2. Charger un fonds de carte par IRIS au format sf
  3. Effectuer une jointure entre les deux fichiers par le code IRIS
  4. Sauvegarder le résultat
  5. Agréger les données statistiques et géométriques par commune
  6. Sauvegarder le résultat.

Préparer les données statistiques

On importe le fichier des individus :

programme

tab_ind<-readRDS("data/menag2018.RDS")

résultat

FALSE   COMMUNE   ARM      IRIS ACHL AEMM
FALSE 1   75056 75101 751010101 C114 2014
FALSE 2   75056 75101 751010101  A11 2008
FALSE 3   75056 75101 751010101  A11 2012

Agréger les données

On commence par créer un tableau long croisant les deux variables et leur effectif pondéré :

programme

tab_long<- tab_ind %>%
           filter(HLML != "Y")%>%
           mutate(COMMUNE = substr(IRIS,1,5)) %>% ### Important ###
    filter(COMMUNE %in% c(75113, 94037, 94041, 94043)) %>%
           group_by(IRIS,HLML)%>%
           summarise(NB=sum(IPONDL))

résultat

IRIS HLML NB
751134901 1 749.88
751134901 2 478.10
751134902 1 763.24
751134902 2 1004.39
751134903 1 81.38

Pivoter le tableau

Puis on fait “pivoter” le tableau pour l’obtenir en format large :

tab_large <- tab_long %>% pivot_wider(id_cols = IRIS, 
                                      names_from = HLML,
                                      names_prefix = "HLM_",
                                      values_from = NB,
                                      values_fill = 0)

résultat

IRIS HLM_1 HLM_2
751134901 749.88 478.10
751134902 763.24 1004.39
751134903 81.38 1118.15
751134904 639.70 221.17
751134905 14.44 1639.81

Ajouter de nouvelles variables

On ajoute de nouvelles variables telles que le nombre total de ménage et le % de ménages en HLM :

tab<- tab_large %>% mutate(TOT = HLM_1+HLM_2,
                           HLM_pct = 100*HLM_1/TOT)

résultat

IRIS HLM_1 HLM_2 TOT HLM_pct
751134901 749.88 478.10 1227.99 61.07
751134902 763.24 1004.39 1767.62 43.18
751134903 81.38 1118.15 1199.52 6.78
751134904 639.70 221.17 860.87 74.31
751134905 14.44 1639.81 1654.25 0.87

Examiner la distribution statistique

On examine l’histogramme donnant distribution statistique du % de ménages ordinaires résidant en HLM par IRIS.

programme

p <- ggplot(tab) + aes (x = HLM_pct) +
                   geom_histogram(breaks = c(0,10,20,30,40,50,
                                             60,70,80,90, 100)) +
                   scale_x_continuous("% de ménages en HLM") +
                   scale_y_continuous("Nombre d'IRIS") +
                   ggtitle(label = "Distribution des logements sociaux dans le 13e arrdt + com. voisines",
                           subtitle = "Source : INSEE, RP 2018")

résultat

Charger les données géométriques

On importe le fichier des iris du Val-de-Marne qui est au format sf en ne gardant que les colonnes utiles

programme

map_iris <- readRDS("data/map_iris.RDS")
map_iris<-map_iris[,c(4,5,1,2,7)] 
names(map_iris)<-c("IRIS","NOM_IRIS","COM","NOM_COM","geometry")
map_iris2 <-map_iris %>%   filter(COM %in% c("75113", "94037", "94041", "94043"))

résultat

FALSE [1] "sf"         "data.frame"
IRIS NOM_IRIS COM NOM_COM
940410101 Plateau Sud 94041 Ivry-sur-Seine
751135007 Gare 7 75113 Paris 13e Arrondissement

Visualisation du fonds iris avec sf

On peut facilement produire une carte vierge des iris du Grand Paris en faisant un plot de la colonne geometry du fichier sf

plot(map_iris2$geometry,col="lightyellow")

Jointure des données IRIS et du fonds de carte

programme

map_iris_tab2<-merge(map_iris2,tab,
                   by.x="IRIS",by.y="IRIS",
                   all.x=T,all.y=F)

résultat

IRIS NOM_IRIS COM NOM_COM HLM_1 HLM_2 TOT HLM_pct geometry
751134901 Salpêtrière 1 75113 Paris 13e Arrondissement 749.88 478.10 1227.99 61.07 MULTIPOLYGON (((652971.2 68…
751134902 Salpêtrière 2 75113 Paris 13e Arrondissement 763.24 1004.39 1767.62 43.18 MULTIPOLYGON (((652861.1 68…
751134903 Salpêtrière 3 75113 Paris 13e Arrondissement 81.38 1118.15 1199.52 6.78 MULTIPOLYGON (((652694 6859…

Sauvegarde du fichier par IRIS

On sauvegarde notre fichier au format .RDS de R

saveRDS(map_iris_tab2,"data/map_iris_hlm2.RDS")

Agrégation statistique + géométriques

Grâce aux nouveaux packages de R (dplyr et sf) il est possible d’agréger simultanément les statistiques et les géométries après les avoir stockés dans un même objet de type “sf”

Du coup, on peut gagner beaucoup de temps dans les traitements et les analyses cartographiques, en particulier si l’on veut tester différents niveaux d’agrégation.

Agrégation des IRIS en communes

L’agrégation est très facile et elle concerne à la fois les variables (de stock) et les geometries

programme

map_com_tab2 <- map_iris_tab2 %>% 
  group_by(COM, NOM_COM) %>% 
  summarise(HLM_1=sum(HLM_1,na.rm=T), 
            HLM_2=sum(HLM_2,na.rm=T)) %>%
  st_cast("MULTIPOLYGON")

map_com_tab2 <- map_com_tab2 %>%  mutate(TOT = HLM_1+HLM_2,
                                  HLM_pct = 100*HLM_1/TOT) 

résultat statistique

COM NOM_COM HLM_1 HLM_2 TOT HLM_pct
75113 Paris 13e Arrondissement 30737 61014 91751 33.5
94037 Gentilly 4117 4716 8833 46.6
94041 Ivry-sur-Seine 9089 19622 28711 31.7
94043 Le Kremlin-Bicêtre 3718 7760 11478 32.4

résultat géométrique

Sauvegarde du fichier par commune

On sauvegarde notre fichier au format .RDS de R

saveRDS(map_com_tab2,"data/map_com_hlm2.RDS")

Visualisation (mapsf)

Le package map_sf

Le package mapsf permet de réaliser des cartes statiques de très haute qualité. Il a en effet été mis au point par des cartographes et des géomaticiens professionnels de l’UMS RIATE. Il prend la suite du package cartography dont la maintenance demeurera assuré quelque temps encore mais ne fera plus l’objet de développements futurs. Le package mapsf présente l’avantage d’être totalement compatibvle avec le package sf ce qui n’était pas autant le cas pour le package cartography, plus ancien, et créé pour être compatible avec l’ancien package sp.

On trouvera la documentation du package mapsf à l’adresse suivante :

https://riatelab.github.io/mapsf/index.html

Création d’un template cartographique

Nous allons dans un premier temps apprendre à créer un fonds de carte vierge mais comportant tout l’habillage nécessaire (“template”). Pour cela nous allons charger différentes couches cartographiques correspondant respectivement au département, aux communes et aux iris :

map_iris<-readRDS("data/map_iris.RDS")
map_com <-readRDS("data/map_com.RDS")
map_dep <-readRDS("data/map_dep.RDS")

map_iris_hlm2<-readRDS("data/map_iris_hlm2.RDS")
map_com_hlm2<-readRDS("data/map_com_hlm2.RDS")

tracé d’un fonds de carte vierge

La fonction mf_map() avec le paramètre type = "base"permet de tracer une carte vide

 mf_map(map_iris_hlm2, type = "base")

Superposition de couches

On peut toutefois ajouter toute une série de paramètres supplémentaire (col=, border=, lwd=, …) et superposer plusieurs fonds de carte avec le paramètre add = TRUE. L’ajout de la fonction layout permet de rajouter un cadre une légende.

# Trace les Iris avec des paramètres
mf_map(map_iris_hlm2,  type = "base", 
       col = "lightyellow", border="gray50",lwd=0.3)
# Ajoute les contours des communes
mf_map(map_com_hlm2,  type = "base", 
       col = NA,border="red",lwd=0.6,
       add = TRUE)

# Ajoute un cadre, un titre et des sources
mf_layout(title = "Paris 13e + voisins", 
          credits = "Sources : IGN et INSEE")

Ajout d’un thème

On peut finalement modifier l’ensemble de la carte en lui ajoutant une instruction mf_theme() qui peut reprendre des styles existants ( “default”, “brutal”, “ink”, “dark”, “agolalight”, “candy”, “darkula”, “iceberg”, “green”, “nevermind”, “jsk”, “barcelona”) mais aussi créer ses propres thèmes

#Choix du thème
mf_theme("darkula")
# Trace les Iris avec des paramètres
mf_map(map_iris_hlm2,  type = "base",
       border="white",
        lwd=0.3)
# Ajoute les contours des communes
mf_map(map_com_hlm2,  type = "base", 
       col = NA, lwd=1,
       add = TRUE)

# Ajoute un cadre, un titre et des sources
mf_layout(title = "Paris 13e & voisins", 
          credits = "Sources : IGN et INSEE")

Ajout de texte

On peut ajouter une couche de texte avec la fonction mf_label(). Par exemple, on va ajouter à la carte précédente le code insee des communes

mf_theme("agolalight")

# Trace les Iris avec des paramètres
mf_map(map_iris_hlm2, 
       type = "base", 
       col = "lightyellow",
       border="gray80",
       lwd=0.3)

# Ajoute les contours des communes
mf_map(map_com_hlm2, 
       type = "base", 
       col = NA,
       border="red",
       lwd=1,
       add = TRUE)

map_iris_hlm2$IRIS_CODE<-substr(map_iris_hlm2$IRIS,6,9)



# Ajoute les codes des communes
mf_label(map_com_hlm2, 
         var="COM",
         cex=1.4, 
         halo = TRUE,
         col="red",
         overlap = FALSE)

# Ajoute les codes des iris
mf_label(map_iris_hlm2, 
         var="IRIS_CODE",
         cex=0.7, 
         col="blue",
         overlap = FALSE)

# Ajoute un cadre, un titre et des sources
mf_layout(title = "Communes et Iris en 2018", 
          frame = TRUE,
          credits = "Sources : IGN et INSEE")

Carte de stock

Une carte de stock représente la localisation de quantités que l’on peut aditionner et dont le total a un sens. Par exemple un nombre d’habitants, un nombre de ménages, un nombre d’automobiles. Ce quantités doivent être représentées par des figures (cercles, carrés, …) dont la surface est proportionelle au stock afin que l’oeil du lecteur puisse les aditionner visuellement.

Dans le package mapsf, on réalise ce type de carte à l’aide de la fonction mf_map()en lui donnant le paramètre type="prop".

On va tenter à titre d’exemple de représenter la distribution du nombre de ménages ordinaires occupant un logement HLM par IRIS :

Carte de stock minimale

Les instructions minimales sont les suivantes :

# Trace les contours des communes
mf_map(x= map_iris_hlm2, 
       type = "base")

# Ajoute le nombre de ménages par IRIS
mf_map(x =map_iris_hlm2, 
      type ="prop",
      var = "HLM_1",
      add=TRUE)

Mais le résultat est peu satisfaisant car les cercles sont trop grands. Il faut en pratique toujours effectuer un réglage de ceux-ci avec l’instruction inches=

Carte de stock habillée

mf_theme("agolalight")
mf_map(map_iris_hlm2, type = "base",  
       col = "lightyellow",border="gray80", lwd=0.3)
mf_map(map_com_hlm2, type = "base", 
       col = NA,border="black",lwd=1,add = TRUE)

mf_map(map_iris_hlm2, var = "HLM_1",type = "prop",
  inches = 0.1, col = "red",leg_pos = "left",  
  leg_title = "Nombre de ménages", add=TRUE)

mf_layout(title = "Distribution des logements HLM en 2018", 
          frame = TRUE,
          credits = "Sources : IGN et INSEE")

Carte choroplèthe

Une carte choroplèthe ou d’intensité représente un phénomène relatif dont la somme n’a pas de sens. Par exemple, il serait absurde d’aditionner les % de logement HLM des IRIS du Val de Marne. Ces variables d’intensité caractèrisent donc l’état général d’une zone (choros) et elles vont être représentées par une couleur appliquée à toute la surface de la zone, d’où leur nom de cartes choroplèthes.

La fonction du package mapsf adaptée aux variables d’intensité est la fonction mf_map()munie du paramètre type = "choro".

On va prendre l’exemple du nombre de voitures par ménage.

Carte choroplèthe minimale

Si on ne précise rien, la carte est réalisée à l’aide de la palette par défaut avec un découpage des classes en quantiles (effectifs égaux).

# Carte choroplèthe
mf_map(
  x = map_iris_hlm2, 
  var = "HLM_pct",
  type = "choro")

Carte choroplèthe habillée

On peut arriver à une carte beaucoup plus satisfaisante en contrôlant l’ensemble des paramètres de couleur et de découpage des classes. Puis en superposant les contours de communes au dessus de la carte des IRIS pour faciliter le repérage.

mybreaks = c(0, 10,20,30,40,50,60,70,80,90, 100)
mypal <- mf_get_pal(n = c(5, 5), pal = c("Greens", "Reds"))
# Carte choroplèthe des iris
mf_map( map_iris_hlm2, var = "HLM_pct",type = "choro",
  breaks = mybreaks,pal = mypal, border=NA,
  col_na = "gray80",leg_title = "% HLM", leg_val_rnd = 0)
# Contour des communes et cadre
mf_map(map_com_hlm2, type = "base", col = NA,
       border="black",lwd=1,add = TRUE)
mf_layout(title = "% de ménages en HLM au RP  2018", frame = TRUE,
          credits = "Sources : IGN et INSEE")

Carte stock + choroplèthe (1)

On peut combiner les deux modes cartographiques par superposition :

mf_theme("agolalight")

# Choisit les classes
mybreaks = c(0,5,10,20,40,80,100)

# Trace la carte choroplèthe
mf_map(
  x = map_iris_hlm2, 
  var = "HLM_pct",
  breaks = mybreaks,
 # pal=mypal,
 type = "choro",
  border="white",
  col_na = "gray80",
 lwd=0.3,
 leg_title = "% ménages", 
 leg_val_rnd = 0,
  
)

# Ajoute les cercles proportionnels

mf_map(
  x =map_iris_hlm2, 
  var = "HLM_1",
  type = "prop",
  inches = 0.06, 
  col = "red",
  leg_pos = "right",  
  leg_title = "Nb ménages",
  add=TRUE
)
# Ajoute les contours des communes
mf_map(map_com_hlm2, 
       type = "base", 
       col = NA,
       border="black",
       lwd=1,
       add = TRUE)

# Ajoute un cadre, un titre et des sources
mf_layout(title = "Les ménages ordinaires en HLM  2018", 
          frame = TRUE,
          credits = "Sources : IGN et INSEE")

Carte stock + choroplèthe (2)

Mais les cercles dissimuent alors les plages de couleur, aussi on peut utiliser le type prop_choro qui place la variable choroplèthe à l’intérieur des cercles

résultat

mf_theme("agolalight")
mybreaks = c(0, 10,20,30,40,50,60,70,80,90, 100)
mypal <- mf_get_pal(n = c(5, 5), pal = c("Greens", "Reds"))
mf_map(map_iris_hlm2, type = "base",  
       col = "gray80",border="white", lwd=0.3)
mf_map(map_com_hlm2, type = "base", 
       col = NA,border="white",lwd=1,add = TRUE)
mf_prop_choro( x = map_iris_hlm2,  var = c("TOT", "HLM_pct"), 
  inches = 0.12, col_na = "grey", pal=mypal,
  breaks = mybreaks, nbreaks = 4, lwd = 0.1,
  leg_pos = c("right", "left"),leg_val_rnd = c(0,0),
  leg_title = c("nb. ménages", "% HLM"),
  add = TRUE)
mf_layout(title = "Les ménages ordinaires en HLM en 2018",
        frame = TRUE, credits = "Sources : IGN et INSEE")

Données spatiales (RPLS, 2020)

Source

Le répertoire des logements locatifs des bailleurs sociaux (RPLS) a pour objectif de dresser l’état global du parc de logements locatifs de ces bailleurs sociaux au 1er janvier d’une année. Il est alimenté par les informations transmises par les bailleurs sociaux. La transmission des informations pour la mise à jour annuelle du répertoire des logements locatifs est obligatoire. Les données sont ensuite géolocalisées à l’adresse et mis à disposition des utilisateurs sur le site du ministère de la transition écologique

Les fichiers sont disponibles en général par régions mais livrés par départements dans le cas de l’Ile de France. Nous allons utilisé ici le fichier du 1er janvier 2020 accessible à l’adresse suivante

https://www.statistiques.developpement-durable.gouv.fr/le-parc-locatif-social-au-1er-janvier-2020-0

Métadonnées

Le fichier de données brutes au format .csv est accompagné d’un document excel précisant le code des variables et la façon dont elles ont été obtenues.

Spatialisation

Le fichier indique pour chaque logement sa localisation précise en terme d’adresse mais aussi d’étage dans un immeuble. A partir de ces données qualitatives, l’INSEE a procédé à un géocodage qui aboutit à la création de deux champs :

  • coordonnées de latitude et longitude non projetées
  • coordonnées de position en projection Lambert officielle

Selon les analyses on peut utiliser l’une ou l’autre de ces coordonnées. Mais la meilleur solution consiste à créer un fichier de type sf (spatial features) en coordonnées WGS94 qu’on pourra ensuite reprojeter dans le système de son choix.

Stratégie

Avant toute exploitation du fichier il est fortement recommandé d’analyser en détail les métadonnées et de définir une stratégie d’analyse.

  1. choisir une première zone d’étude de petite taille et localisée de préférence dans un espace que l’on connaît bien.
  2. choisir des variables intéressantes dont l’on connaît bien la signification et dont on a analysé en détail les métadonnées
  3. vérifier la qualité des données en regardant notamment le nombre de valeurs manquantes, le dégré de précision, etc.
  4. sélectionner des données auxiliaires issues d’autres sources que l’on souhaite croiser avec celles du RPLS en s’assurant de leur compatibilité (espace, temps, définition, …)
  5. Ajouter les coordonnées spatiales et stocker le résultat dans un fichier de type sf comportant les indications de projection.

Importation du fichier

On importe le fichier enregistré au format RDS et on vérifie sa taille avec dim() et sont ype avec class()

don <- readRDS("data/RPLS2020.RDS")
dim(don)
[1] 844302     73

Le tableau comporte 844302 lignes (chacune correspondant à un logement) et 73 variables (décrites dans les métadonnées).

Choix de la zone d’étude

On décide de limiter notre analyse dans un premier temps au 13e arrondissement et aux trois communes voisines.

Choix des variables

On va se limiter ici à un très petit nombre de variables

variables de localisation

  • result_id : code de l’adresse
  • result_label : label de l’adresse
  • LIBCOM : nom de la commune
  • DEPCOM : code de la commune
  • latitude : coordonnées latitude
  • longitude: coordonnée longitude
  • X : coordonnée projetée (EPSG = 2154)
  • Y : coordonnée projetée (EPSG = 2154)

variables thématiques

  • CONSTRUCT : année de construction
  • SURFHAB : surface habitable en m2
  • NBPIECE : nombre de pièces

Extraction du fichier

On applique la double sélection des individus et des variables en nous servant des fonctions filter()et select()du package dplyr.on aboutit ici à un fichier de 8139 lignes et 11 variables.

sel <- don %>% mutate(DEPCOM = substr(result_id,1,5)) %>% ### Important ###
  filter(DEPCOM %in% c("75113", "94037", "94041", "94043"))%>%
  select(result_id, result_label,
         DEPCOM, LIBCOM, 
        latitude,longitude,X,Y,
        CONSTRUCT,SURFHAB,NBPIECE) 

dim(sel)
[1] 54419    11

Recodage et typage

Certaines variables doivent être recodées ou changées de type afin de faciliter leur exploitation ultérieure par R.

sel$DEPCOM <- as.character(sel$DEPCOM)
sel$LIBCOM <- as.factor(sel$LIBCOM)
sel$PLG_IRIS <- paste(sel$DEPCOM,sel$PLG_IRIS, sep = "")
sel$SURFHAB <- as.numeric(sel$SURFHAB)

Résumé rapide

On analyse rapidement les variables thématiques choisies

CONSTRUCT SURFHAB NBPIECE
Min. :1850 Min. : 8.00 Min. :1.000
1st Qu.:1958 1st Qu.: 42.00 1st Qu.:2.000
Median :1971 Median : 57.00 Median :3.000
Mean :1971 Mean : 56.87 Mean :2.744
3rd Qu.:1990 3rd Qu.: 72.00 3rd Qu.:4.000
Max. :2019 Max. :200.00 Max. :9.000

Sauvegarde du fichier

On sauvegarde le fichier obtenu au format .RDS afin de garder le formatage des variables :

saveRDS(sel,"data/sel_logt2.RDS")

Localisation spatiale

Retour sur sf

Nous revenons sur le package sf (spatial features) que nous avons déjà rencontré au moment de la création de cartes thématiques par IRIS ou communes à l’aide du package mapsf.

Ici le package sf va être utilisé pour cartographier non pas des zones mais des localisations ponctuelles. Il pourra être à nouveau couplé avec le logiciel de cartogaphie statique comme mapsf , afin par exemple de placer les localisations des logements sociaux au dessus du fonds de carte des IRIS ou communes.

Mais il pourra aussi servir de base à des cartographies dynamiques permettant de placer les points sur des réseaux de rue et plus généralement sur des “tuiles” cartographiques permettant d’effectur des zoom. On utilisera à cet effet d’autres packages comme leaflet ou sa version simplifiée mapview.

Données ponctuelles

Nous reprenons le fichier de localisation établi au chapitre précédent et nous ne conservons que 6 variables:

logt <- readRDS("data/sel_logt2.RDS") %>%
        select(adresse=result_id,
               X,Y,
               date = CONSTRUCT)
adresse X Y date
75113_6693_00078 653515.8 6858687 1987
75113_6693_00078 653515.8 6858687 1987
75113_6693_00078 653515.8 6858687 1987

Données IRIS

Nous chargeons par ailleurs le fichier des IRIS relatif à la zone d’étude :

map_iris <- readRDS("data/map_iris_hlm2.RDS") 
IRIS NOM_IRIS COM NOM_COM HLM_1 HLM_2 TOT HLM_pct
751134901 Salpêtrière 1 75113 Paris 13e Arrondissement 749.88444 478.1026 1227.987 61.066153
751134902 Salpêtrière 2 75113 Paris 13e Arrondissement 763.23525 1004.3883 1767.624 43.178609
751134903 Salpêtrière 3 75113 Paris 13e Arrondissement 81.37525 1118.1495 1199.525 6.783958

Agrégation par commune

Rappel : on peut agréger les géométries d’un fonds sf. Ici on va créer le fonds de carte des communes.

map_com <- map_iris %>% group_by(COM,NOM_COM) %>%
                summarise() %>%
                st_cast("MULTIPOLYGON")

Vérification de la projection

Nous savons que les coordonnées X,Y du fichier logement sont projetées en EPS 2154. Mais quelle est la projection de notre fonds IRIS ? S’agit-il de la même ?

st_crs(map_iris)$proj4string
[1] "+proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
st_crs(2154)$proj4string
[1] "+proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

A priori il s’agit bien de la même de sorte que les coordonnées X,Y devraient bien se superposer sur le fonds IRIS

Test de superposition

par(mar=c(0,0,0,0))
#trace les iris
plot(map_iris$geometry, 
     col="lightyellow", border="gray70",
     lwd=0.2)
# trace les communes     
plot(map_com$geometry, 
     col=NA, lwd=1, add=T)
# ajoute les points
points(x=logt$X,
       y=logt$Y, 
       cex=0.2,
       col="red",
       pch = 16)

fichier des adresses

Nous allons maintenant établir un fichier de localisation des adresses en nous servant de l’identifiant unique fourni par l’INSEE.

adr <- logt %>% select(adresse,X,Y) %>% 
               filter(duplicated(adresse) == F) %>%
               filter(is.na(X) ==F,is.na(Y)==F)

On constate qu’il n’y a que 652 adresses différentes alors que notre fichier fait état de 8139 logements. Une adresse regroupe donc en moyenne plus de 10 logements (habitat collectif).

Transformation en fichier sf

La transformation de notre fichier initial au format sf est facile à réaliser avec la fonction st_as_sf() du package sf. Mais il faut prendre garde de bien préciser le système de projection si l’on veut pouvoir ensuite l’utiliser.

map_adr <- st_as_sf(adr, coords = c("X","Y"))
st_crs(map_adr)<- 2154
str(map_adr)
Classes 'sf', 'data.table' and 'data.frame':    1729 obs. of  2 variables:
 $ adresse : chr  "75113_6693_00078" "75113_1918_00006" "75113_2852_00001" "75113_6589_00006" ...
 $ geometry:sfc_POINT of length 1729; first list element:  'XY' num  653516 6858687
 - attr(*, "sf_column")= chr "geometry"
 - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA
  ..- attr(*, "names")= chr "adresse"

Agrégation des logements

Notre nouveau fichier sf permet désormais d’effectuer des jointures avec le fichier des logements sociaux. A titre d’exemple on peut désormais compter le nombre de logements par adresse et leur ancienneté moyenne.

programme

logt_by_adr <- logt %>% 
               group_by(adresse) %>%
               summarise(nblog = n(),
                         datemoy = mean(date))

résultat

adresse nblog datemoy
75113_0027_00003 1 1996
75113_0027_00005 1 1996
75113_0027_00007 1 1996
75113_0027_00011 34 1997
75113_0027_00013 19 1997
75113_0027_00015 20 1997
75113_0027_00017 34 1997
75113_0028_00025 12 1965
75113_0028_00027 13 1965
75113_0028_00029 13 1965

Jointure

On peut désormais effectuer la jointure entre les données agrégées par adresse et le fichier sf de localisation des adresses :

map_logt <- inner_join(logt_by_adr,map_adr) %>% st_as_sf()

Cartographie avec mapsf

On peut désormais utiliser les méthodes de cartographie déjà vues avec mapsf :

mf_theme("agolalight")
mybreaks = c(1900, 1950, 1960, 1970, 1980, 1990, 2000, 2010, 2020)
library(RColorBrewer)
mypal=brewer.pal(n = 8,name = "Spectral")
mf_map(map_iris, type = "base",  
       col = "gray80",border="white", lwd=0.3)
mf_map(map_com, type = "base", 
       col = NA,border="black",lwd=1,add = TRUE)
mf_prop_choro( x = map_logt,  var = c("nblog", "datemoy"), 
  inches = 0.1, col_na = "grey", pal=mypal,
  breaks = mybreaks, nbreaks = 4, lwd = 0.1,
  leg_pos = c("right", "left"),leg_val_rnd = c(0,0),
  leg_title = c("nb. logements", "ancienneté"),
  add = TRUE)
mf_layout(title = "Les logements sociaux en 2020",
        frame = TRUE, credits = "Sources : IGN et RPLS")

Sauvegarde des fichiers cartographiques

On sauvegarde nos différents fichiers cartographiques au format sf relatifs à la zone d’étude.

saveRDS(map_com,"data/sel_map_com2.RDS")
saveRDS(map_iris,"data/sel_map_iris2.RDS")
saveRDS(map_logt,"data/sel_map_logt2.RDS")

Cartographie dynamique

Statique ou dynamique ?

  • Cartographie statique
    • production d’images fixes de qualité
    • respect strict des règles de la sémiologie graphique
    • choix libre d’une projection adaptée (e.g. EPSG 2154)
    • production de documents imprimés à finalité normative ou scientifiques
  • Cartographie dynamique
    • production d’interfaces consultables dans un navigateur.
    • modification possible de l’échelle et de l’arrière-plan
    • projection imposée par les “tuiles” (EPSG 4326)
    • production de documents interactifs à finalité citoyenne ou exploratoire

Packages R de cartographie dynamique

  • leaflet : la référence
    • Une librairie javascript non liée à un langage (R, Python, html, …)
    • Disponible dans R sous forme de package
    • Développement constant
  • ggmap : l’empire contre attaque
    • des outils cartogaphiques utilisant la syntaxe de tidyverse
    • impose désormais un lien avec Google
  • tmap : une solution hybride
    • permet de passer facilement du mode statique au mode dynamique
  • mapview : l’équivalent de mapsf
    • mis au point par des développeurs allemands
    • facilite l’usage de leaflet
    • en progrès constant (mais instable)

Préparation des données

On charge les fichiers au format sf et on les transforme en projection WGS94 (EPSG=4326), condition indispensable pour ajouter des “tuiles” dynamiques lors des zoom.

map_com <- readRDS("data/sel_map_com2.RDS") %>%
              st_transform(4326)
map_iris <- readRDS("data/sel_map_iris2.RDS") %>%
              st_transform(4326)
map_logt <- readRDS("data/sel_map_logt2.RDS") %>%
              st_transform(4326)

Carte par défaut

Mapview produit par défaut une carte dynamique du fichier sf.

mapview(map_logt)

Superposition de couches

On peut créer des couches et les aditionner avec ‘+’ :

m1 = mapview(map_com, zcol = "NOM_COM") 
m2 = mapview(map_logt)
m1+m2

Exemple complet

On va essayer de reproduire la carte statique faite avec mapsf

# Carte des communes
map1 <- mapview(map_com, lwd=1, legend= FALSE,
                alpha.regions = 0.1)
# Carte des iris
map2 <- mapview(map_iris,lwd = 0.3, label= "NOM_IRIS",
                legend= FALSE, alpha.regions = 0)
# Carte des logements
map3 <- mapview(map_logt,
                zcol = "datemoy",
                at = c(1900,1960, 1970,1980,
                       1990,2000,2010, 2021),
                col.regions = brewer.pal(8, "Spectral"),
                cex= "nblog")
map1+map2+map3